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. Abstract. Two-dimensional simulations of strongly anisotropic supernova explosions of a nonrotating 15 M© blue supergiant 

■ progenitor are presented, which follow the hydrodynamic evolution from times shortly after shock formation until hours later. 
It is shown that explosions which around the time of shock revival are dominated by low-order unstable modes (i.e. by a 

■ superposition of the 1-2 and I - \ modes, in which the former is strongest), are consistent with all major observational 
^ ' features of SN 1987 A, in contrast to models which show high-order mode perturbations only and were published in earlier 

, work. Among other items, the low-mode models exhibit final iron-group velocities of up to ~ 3300 km/s, strong mixing at 

the He/H composition interface, with hydrogen being mixed downward in velocity space to only 500 km/s, and a final prolate 
anisotropy of the inner ejecta with a major to minor axis ratio of about 1.6. The success of low-mode explosions with an energy 
of about 2 X 10^^ erg to reproduce these observed features is based on two eff'ects: the (by 40%) larger initial maximum velocities 
of metal-rich clumps compared to our high-mode models, and the initial global deformation of the shock. The first eff'ect protects 
\ the (fastest) clumps from interacting with the strong reverse shock that forms below the He/H composition interface, by keeping 

their propagation timescale through the He-core shorter than the reverse shock formation time. This ensures that the outward 
motion of the clumps remains always subsonic, and that thus their energy dissipation is minimal (in contrast to the supersonic 
case). The second eff'ect is responsible for the strong inward mixing of hydrogen: The aspherical shock deposits large amounts 
Q ■ of vorticity into the He/H interface layer at early times (around t - 100 s). This triggers the growth of a strong Richtmyer- 

%^ I Meshkov instability that results in a global anisotropy of the inner ejecta at late times (i.e. around t - 10000 s), although the 

C/) , shock itself has long become spherical by then. The simulations suggest a coherent picture, which explains the observational 

data of SN 1987 A within the framework of the neutrino-driven explosion mechanism using a minimal set of assumptions. It 
is therefore argued that other paradigms, which are based on (more) controversial physics, may not be required to explain this 
event. 
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1. Introduction 



Key words, hydrodynamics - instabilities - nucleosynthesis - observed features in a satisfactory way. An extensive bibliog- 
shock waves - stars: supernovae raphy and overview of such work was given in the first paper 

of the present series (Kifonidis et al. 2003, henceforth Paper I). 

In this latter reference we ourselves have embarked on a 
somewhat more ambitious effort, namely to investigate the ob- 
Large-scale anisotropics and mixing, as well as smaller-scale servational consequences of the modem (i.e. multidimensional, 
clumping, are common features of supernova (SN) remnants o^" convectively supported) neutrino-driven explosion paradigm 
like Cas A, or SN 1987 A, and have been attributed to the occur- for ^ore collapse supernovae. For this purpose we performed 
rence of hydrodynamic instabilities (in particular the Rayleigh- ^n exploratory, high-resolution 2D simulation of a type II su- 
Taylor or RT instability) in the SN explosions that gave rise P^mova in a 15 Mo blue supergiant which included a detailed 
to these remnants. Starting the supernova by depositing kinetic modellmg of the explosion itself, and which we dubbed "Model 
and/or thermal energy into a progenitor model, many groups T310a". This calculation covered the time from 20 ms up to 
in the past have performed multidimensional hydrodynamic ^ore than 5 hours after core bounce, and, at the time of sub- 
simulations to study the mixing associated with these effects, ^ssio" of *e present paper, is still the only multidimensional 
However, so far the models have not been able to explain the simulation that followed the hydrodynamic evolution and the 

growth of anisotropics from the earliest moments after core 

Send offprint requests to: K. Kifonidis collapse until well after the time of shock emergence from the 

Correspondence to: kok@mpa-garching.mpg.de stellar photosphere. In this way the effects of neutrino heating, 
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convection, and explosive nucleosynthesis deep inside the core, 
as well as their impact on the Rayleigh-Taylor instabilities in 
the stellar envelope were all taken into account. 

For simulating the onset of the explosion in Paper I, we 
approximated the neutrino transport by a simple light-bulb 
scheme, and parametrized the neutrino fluxes emitted by the 
optically thick proto-neutron star core. In an attempt to repro- 
duce the high nickel velocities of SN 1987 A, the most sensi- 
tive test case for supernova theory available to date, we per- 
formed the calculation using rather high initial values for the 
neutrino fluxes, and a subsequent exponential decay of these 
fluxes with time. This prescription initially appeared promis- 
ing, since it gave a rather energetic explosion. However, the 
maximum nickel velocities that were finally obtained from the 
simulation, were disappointingly small. 

Furthermore, because the explosion in this model set in 
very rapidly, it did not allow hydrodynamic instabilities, which 
occur in the neutrino-heated post-shock layers, to grow to suf- 
ficient strength such that they could deform the shock. Hence 
the model showed only relatively small-scale deviations from 
spherical symmetry - i.e., higher order modes of the fluid flow 
in terms of an expansion of the inhomogeneities in Legendre 
polynomials of order / - and was not able to reproduce the 
global anisotropy of the ejecta and the extent of the mixing 
observed in SN 1987 A. It also left behind a small neutron star 
of only about 1.1 M© (baryonic). Apparently, our parametriza- 
tion of the neutrino fluxes had resulted in too short an explosion 
timescale, and our hope was that the aforementioned problems 
might be overcome with an improved description of the neu- 
trino efl'ects, e.g. by a modification of the (ad hoc) assumptions 
on the behavior of the boundary conditions for the neutrino pa- 
rameters. 

This provided the motivation for the systematic two- 
dimensional neutrinoji^rod^ study that we 
presented in IScheck et aiJ ( l2nn4l2nn4 . In these works, we re- 
placed the high initial values of the core neutrino fluxes with 
lower ones, and the (old) exponential decay with a much slower 
decline over time, an assumption which is in agreement with 
the results of current simulations employing Boltzmann neu- 
trino transport. (Note that the core luminosity is imposed at our 
inner grid boundary which follows a Lagrangian mass coordi- 
nate.) Furthermore, we dropped the simple neutrino light bulb 
approximation and replaced it with a new, grey, characteristics- 
based neutrino transport scheme that we use in the free- 
streaming and semi-transparent regimes. This allows us to take 
into account the efl'ects of neutrino heating and cooling on the 
flux, e.g. to include the efl'e ct of mass accretion on the neutrino 
luminosities (for details see Scheck et al. 2006). 

With this new approach Scheck et al. (20041 12004 were 
able to demonstrate that, if the explosion sets in sufficiently 
slowly, low-mode hydrodynamic instabilities, i.e. low-mode 
convection ( Chandrasekhai 1961: Her ant 1995) in the neutrino- 
heating layer, and the recently discovered advective-acoustic 
(FoeHzzo & Taeeer "2000"; FoeHzzo "2002"; "Foglizzo & Galettil 
.2003; Fog lizzo et al. 2005; Ohriishi et al. 2005), and SASI in - 
stabilities (|Blondin et alJ|2003|;|Blondin & Mezzacappa|2005|), 
can grow from small random perturbations behind the shock 
and result in a global anisotropy of the shock and the ejecta. 



Moreover such explosions with typical supernova energies also 
leave behind neutron stars with reasona ble masses a nd h igh re- 
coil velocities, all at the same time. Sche ck et"all |2006) fur- 
thermore showed that the occurrence of different low-order un- 
stable modes (/ = 1, / = 2) in the post- shock ffow, might even 
explain a possible bimodality of the observed galactic pulsar 
velocity distribution without the need to assume any physics in 
addition to what is already part of the neutrino-driven explosion 
paradigm. 

Yet, it needs to be pointed out that currently the viability 
of neutrino-driven supernova explosions for progenitor stars 
of more than ~ IOM^t) is still an open question (see, e.g., 
iBuras et aljEooi l2006alhh. Our adherence to the neutrino- 
driven explosion mechanism and our procedure of triggering 
supernova explosions by the use of a boundary condition for 
the neutrino ffux may be justified by the fact that uncertain- 
ties are still present even in the most sophisticated recent 2D 
supernova simulations, despite the significant improvements in 
the treatment of neutrino physics and transport that have been 
achieved compared to the first generation of multidimensional 
modelling. Such recent 2D simulations, which include spectral 
neutrino transport, (see Buras et aL .2003.. 2006ah) could not 
confirm the success of first generation models, whi ch employed 
gray neutrino diffusion and found explosions (e.g. Herant et alJ 
j1 994i_Burrows et aLJ 995: .Frver.J999: Jrver & Warren 20021 

However, it should be noted that 3D simulations with spec- 
tral transport have not been performed yet, and might reveal 
important differences in the flow dynamics compared to ax- 
isymmetric models, in particular with respect to the growth of 
convective and Rayleigh-Taylor instabilities, non-radial accre- 
tion shock instabilities, and vortex behaviour. Moreover, the 
question whether convective or mixing processes below and 
around the neutrino sphere could enhance the neutrino spheric 
neutrino emission and thus support the neutrino heating behind 
the shock, must still be considered as unsettled: Doubly dif- 
fusive i nstabilities (Bruenn et al. 2006), neutrino-bubble insta- 
biliti es Js ocrates e t alJl2005l) . or magnetic buoyancy instabili- 
ties (Wilson et al. z005) deserve further investigation by mul- 
tidimensional modelling. Also magneto-rotational instabilities 
were suggested to aid the explosion by creating viscous heating 
behind the shock in addition to the still dominant energy input 
there by neutrinos (Thom pson et al. 2005) . 

Although these currently unre solved issues need to be kept 
in mind, the (2D) simulations of "Scheck et al." ("2006") suggest 
a quite remarkable perspective, namely that many fundamen- 
tal properties of observed supernovae and neutron stars might 
be traced back to the same origin, i.e. to non-radial hydrody- 
namic instabilities during the first second of a neutrino-driven 
SN explosion. However, as these calculations did not cover also 
the later phases of the explosion, they could not provide in- 
sight into the inevitable interaction of the inner ejecta with the 
(possibly massive, e.g. hydrogen) envelope of the progenitor. 
Since this interaction (which happens on a timescale of hours 
to days, depending on the type of the progenitor) determines 
many of the more intricate observational features, especially of 
Type II supernovae, an important question still remains unan- 
swered: Can globally anisotropic neutrino-driven explosions 
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also lead to high final nickel velocities, strong mixing at the 
He/H composition interface, large final ejecta anisotropics, and 
hence sizeable polarization of the emitted light? 

A positive answer to this question would not only off'er a 
(natural) explanation for the gamut of observational data col- 
lected in case of SN 1987 A, or SN remnants like Cas A. It 
would also provide substantial support to the neutrino mech- 
anism. The present paper is an attempt to explore this ques- 
tion in the framework of two-dimensional, i.e. axisymmet- 
ric simulations. We exemplarily follow the evolution of glob- 
ally anisotro pic explosions of the type studied in Sc heck et alJ 
(l2004 E 006) from several milliseconds up to several hours af- 
ter core bounce in high-resolution simulations. 

Although we do not simulate the "full SN problem", but 
trigger the neutrino-driven explosion by the use of a suit- 
able inner boundary condition, we nevertheless think that our 
modelling approach is adequate to address the aforementioned 
points because it is "consistent" in the following sense: The dis- 
tribution of kinetic and internal energy in the flow, the number 
and size of flow structures due to hydrodynamic instabilities, 
the shape and asymmetry of the supernova shock, and the dis- 
tribution of nucleosynthetic products (both in real and velocity 
space) develop naturally in a detailed neutrino-hydrodynamic 
calculation, in response to a boundary condition for the neu- 
trino flux and the contraction of the neutron star core. Neither 
of the former quantities, which determine the kind of observ- 
ables that we are interested in, is put in "by hand" - in contrast 
to previous work in which ad hoc prescriptions were used for 
at least one of these aspects. 

The question of how close our boundary condition and initi- 
ation of the explosion come to reality cannot be finally assessed 
at present, because of the lack of detailed observed neutrino 
light curves, and/or the aforementioned uncertainties inherent 
to more sophisticated modelling of the (neutrino-driven) explo- 
sion mechanism. Yet, the choice of particular parameter values 
for this boundary condition aff'ects mainly the integral charac- 
teristics of the explosion like its total energy and time of onset, 
as well as the integral properties of the n eutron star, e.g. its 
mass and contraction time scale (see again Sch eck et al.ll200(j 
for details). A particular parametrization is in the end vindi- 
cated by its success to reproduce the observational features of 
SN 1987A. 

The stellar model that we shall c onsider is the ISM© 
blue supergiant SN 1987 A progenitor of lWooslev et al.l (Il988h . 
which we already used in Paper I. This makes the present work 
a direct sequel and extension of this earlier investigation. It is 
our hope that these papers, together with the constraints ob- 
tained from observations, will serve to elucidate the still poorly 
understood hydrodynamics of SN 1987 A, and help to disclose 
the true nature of this and similar events. 

We shall proceed with a description of the most recent mod- 
ifications to our computer codes in Sect. |2l We then comment 
briefly on our initial data and some computational aspects in 
Sect. 13 The results of our calculations are presented in Sect.lU 
and compared with SN 1987 A in Sect. 13 Our conclusions can 
be found in Sect. |6l The appendix finally contains a su mmary 
of the AUSM+ numerical flux function of lLioul (Il996l) that is 
used in one of our hydrodynamics codes. 



2. Modifications to the codes 

The computer programs that were used for the present calcu- 
lations , are the neutrino hydrodynamics code of IScheck et al.l 
2006) and the version of the adaptive mesh refine- 
ment hydrodynamics code Amra that we described in Paper I. 
Both use the direct Eulerian versi on of the Piec ewise Parabolic 
Method (PPM) of C olella & Woodward! (Il 984 to integrate the 
hydrodynamics equations. Several changes were made to the 
former code, though. In particular we have improved its poten- 
tial to resolve shocks, by updatin g the hybrid R i emann solve r 
algorithm which we employed in IScheck et al.l (l2004 l2006h . 
In these works we have made use of the standard (in PPM) ex- 
act Riemann solver to evolve all zones on the grid, except for 
those zones which are located in the vicinity of strong (grid- 
afigned) shocks, for which the HLLE flux of Einfeldt ( 198^ 
was employed (cf . Paper I). A fter replacing the HLLE flux by 
the AUSM+ flux of Lioul(ll99 6) (see also AppendixlAl we now 
achieve both, post-shock states completely free of odd-even de- 
coupling, and very sharp shocks with only one or two interior 
points (while with the HLLE flux at least 3-4 points are re- 
quired). 

We have also completely revised the nucleosynthesis mod- 
ules of this code (see the end of this section for our motiva- 
tion). For zones with temperatures 1.5xlO^K<r< T^^^ 
we now solve a 32 species nuclear reaction network instead 
of the modified a-nuclei network described in Paper I. The iso- 
topes considered in the 32 species network are listed in TableQ 
Nuclear data and rates for 141 reactio ns in total stem from the 
compilation of F.-K. Thielemann (see Thielemann et al]ll996t 
Hoff'manetal. 1999, and the references therein). The ordi- 
nary diff'erential equations (ODEs) constituting the nuclear net- 
work are solved with the variable-order Bader-Deuflhard semi- 
implicit, stifl' ODE integrator, which employs sub-stepping 
and Richardson extrapolation to achieve very large eff'ec- 
tive (time) ste ps (|Bader & Deuflhard|[T983[ [Press et aLlEHl 
|Timme|l999|). 

For T > T^^^ a nuclear statistical equilibrium (NSE) solu- 
tion is computed: Given the temperature, T, density, p, electron 
fraction, Yq, and the parameters Z/ (charge number). A/ (mass 
number), Bi (binding energy) and Gi(T) (partition function) of 
species /, we iteratively solve (with a globally convergent mul- 
tidimensional Newton-Raphson scheme) the Saha equations 

(1) 

and the equations of charge and mass conservation 

J^ZiYi = Ye, YjAiYi = l (2) 

for the neutron and proton abundances, and Fp, respectively. 
From these follow the abundances 7/ of all other species by 
virtue of Eq. ([l}. 

Nucle ar screening was included in the NSE solution ac- 
cording to lHix & Thielemannl(ll996l) . but is currently neglected 
in the nuclear network solver. In order to get a perfect match 
between the numerical solution of the nuclear network ODEs 
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Table 1. Isotopes of the 32 species nuclear reaction network 



n 


P 


^He 




160 


20Ne 




2«Si 




36 


«Ca 


44 








55Cr 




55 Mn 


5*Mn 


«Fe 




55Fe 


S'^Fe 


"Fe 


58Fe 


55Co 




56Ni 


"Ni 


58Ni 


60Ni 





and the NSE solution we had to set = 8 x 10^ K, i.e. 
we use the nuclear network well into the medium-temperature 
NSE regime, where the composition consists of neutrons, pro- 
tons, and a-particles, but essentially no heavy nuclei. The re- 
sulting algorithm proved to be very robust and shows no dif- 
ficulties in making smooth transitions from temperatures as 
high as ~ 10^^ K to less than 10^ K, and vice versa. In ad- 
dition, it turned out to be very efficient. Only a small number 
of (main) steps per burning zone are necessary with the Bader- 
Deuflhard method to obtain a relative integration accuracy of 
10"^. The employed T^^^ turned out to be still low enough to 
avoid the reaction rates becoming extremely large, and hence 
nuclear timescales becoming extremely small ^ 

The coupling of the nuclear burning algorithm with the 
hydrodynamics was done in the same operator- split, "co- 
processing" approximation that we described in Paper I, i.e. 
energy and composition changes from the 32 species solution 
are currently not fed back into the solution of the hydrody- 
namics, which instead still relies on the smaller 4-species NSE 
solution tabulated in the equation of state of Janka & Miiller 
( 1996). Work is in progress to extend the implicit nuclear algo- 
rithm with an energy equation (see Miiller 1986), and with the 
leptonic equation of state of Timmes & Swesty (2000), in order 
to obtain a fully coupled treatment of the composition changes, 
nuclear source terms, equation of state, and hydrodynamics in 
future simulations. 

We note here that our primary motivation for coding the 
present method was its enhanced robustness and efficiency 
over the implementation of the nuclear physics that we used 
in Paper I. This enables us to handle situations in which the 
previous algorithm (wich was based on the "classical" first- 
order Euler backward time discretization), either experienced 
convergence problems or required too many nuclear time steps 
to cover a hydrodynamic step. As a side eff'ect, and due to the 
large gain in efficiency provided by the Bader-Deuflhard inte- 
grator (see iTimmesll 19991) . we can now also solve larger reac- 
tion networks. However, improvements in the accuracy of the 
nuclear yields which may thereby be achieved, need to be con- 
trasted with the fact that the relative yields of many isotopes 
within the iron group (e.g. the ^^Ni/^^Fe ratio), depend ex- 
tremely sensitively on the electron fraction, Yq, and thus require 
also a more sophisticated, non-grey neutrino transport scheme 
for their accurate computation than the neutrino treatment we 
employ here (see also the discussion of these issues in Paper I). 
For this reason we will give only integrated abundances (or 

^ Despite our efficient solution of the 32-species network, the sim- 
ulations still require about a factor of five more computer time than 
without computing the nucleosynthesis. 



mass fractions) for the iron group in this paper, instead of indi- 
vidual abundances for single iron-group nuclei. 

3. Computational aspects 

3.1. Numerical mesh, initial, and boundary conditions 

Throughout this work, 2D spherical coordinates (r, ?^), a com- 
putational grid with 180° width in and axisymmetry are 
adopted. To perform the present simulations we chose to stick 
to the same ISM© blue supergiant model of IWooslev et alJ 
('l988) and the post-bounce model of'Bruenn (1993*) (at a time 
of 20 ms after core bounce) that we have already used as initial 
data in Paper I. 

Similar to this latter work, we follow the first second of the 
explosion using neutrin o-hydrodynamics ca lculations. These 
are set up as described in lScheck et alJ (l2006h and we refer the 
reader to this paper for details, in particular for the employed 
boundary conditions. 

We wish to highlight one point though, which is connected 
to the fact that we start our 2D simulations from supernova 
progenitor and post-bounce models that are one-dimensional. 
Unlike the grid-less smooth particle hydrodynamics (SPH) 
method, our grid-based Eulerian hydrodynamics codes exhibit 
such a low level of numerical noise that a one-dimensional, 
isotropic initial configuration remains perfectly isotropic when 
evolved on a 2D grid, even if, e.g., a stratification is present 
in the ID initial data which is prone to become convective. To 
trigger the growth of non-radial hydrodynamic instabilities in 
the post-shock flow, we thus need to explicitly apply a small 
random perturbation to the initial model, which we add to the 
velocity field and for which we typically use an amplitude of 
0.1%. This perturbation is applied only once, at the beginning 
of the calculation, i.e. at ^ 20 ms post-bounce. 

To prescribe the neutrino radiation field of the contracting 
neutron star core we make use of the same inner boundary con- 
traction and parametrization (with A^^^^.^ = 0.1 8MnC^,r^, = 
1.0 s) that was employed for Model B18 in IScheck et alJ 
(l2006h . This resulted in a total electron flavour neutrino lumi- 
nosity of the neutron star core, L^^^^ = Ly^ + Ly^ = 4.45 x 
10^^ erg/s, which was assumed to be constant during the first 
second. 

In order to explore diff'erent explosion asymmetries result- 
ing from the non-radial instabilities, we calculated two models 
with this particular setup that diff'er only in the initial pertur- 
bation. To achieve this we simply used two diff'erent random 
number sequences to perturb the initial velocity field in these 
two runs. Everything else, in particular the perturbation ampli- 
tude that we employed, was kept the same. The resulting mod- 
els will henceforth be called bl8a and bl 8b. We point out here 

- and elaborate on this in more detail in ' Scheck et alJ (l2006h 

- that already such minor diff'erences between the simulations 
produce diff'erent ejecta asymmetries, because the growth of the 
initial perturbations is extremely nonlinear. In fact it is chaotic. 

We also used a second parametrization for the inner bound- 
ary condition which diff'ers from the former setup only in a 
somewhat larger gravitational binding energy loss of the neu- 
tron star core, A^^^e = 0.23 MqC^. This leads to an L^°^^ which 
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entropy [kg/nuc] entropy [kg/nuc] entropy [kg/nuc] 

5 10 15 20 25 30 5 10 15 20 25 5 10 15 20 25 30 




0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 

X [10'^ cm] X [10'^ cm] x [10^^ cm] 



Fig. 1. Top row: Entropy distribution (in units of ^e/nucleon) for three of the simulations listed in Table |2| at early times. From 
left to right a) Model bl8a at ^ = 1 s, b) Model bl8b at ^ = 1 s, and c) Model b23a at t = 0.92 s. Bottom row: Density distribution 
for the same models at later instants. From left to right: d) Model bl8a at ^ = 3024 s, e) Model blSb at ^ = 3031 s, and f) Model 
b23a at ^ = 3083 s. Note the change of shape of the supernova shock (the outermost discontinuity in all plots). Note also the polar 
jets which deform the otherwise spherical main shock at late times. They are a consequence of our use of a spherical coordinate 
grid and our assumption of axisymmetry. 
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is by 28% higher than in the bl8 models, and gives rise to 
Model b23a, which attained an explosion energy that is twice 
as large as that of the bl8 case (see Table |3 and Sect. 14. 111 . 

Note that throughout this paper we use lower-case initial 
letters in our model names to distin guish the present simula- 
tions from the models of iScheck et al.l (l2006l) . which did not 
make use of the nucleosynthesis solver, and were still per- 
formed using the HLLE flux instead of the AUSM+ flux (see 
above). Note also that in all of our simulations (which were 
carried out with 400 radial and 180 angular zones), the inner 
boundary condition for neutrinos was assumed to be isotropic, 
and that the core luminosities were held constant over the 
time ti (for a motivation and full details, including a compar- 
ison with results f rom Boltzmann transport calculations, see 
[SchecketalJl200(J). 

The hydrodynamic evolution after the first second (which 
includes the growth of Rayleigh-Taylor and Richtmyer- 
Meshkov instabilities) was simulated with the Amra code 
and the mesh refinement hierarchy, mesh resolution, remap- 
ping scheme, and boundary conditions that were described in 
Paper I. This allows us to temporarily achieve a resolution 
equivalent to covering the entire progenitor with 2.6 x 10^ 
equidistant radial zones, although the eff'ective number of ra- 
dial zones on the finest level of the AMR grid hierarchy was 
only 3072. The eff'ective number of angular zones on this re- 
finement level was 768. 

3.2. Constraints due to the grid-geometry 

At the outset we must acknowledge that a two-dimensional 
hydrodynamic study in standard spherical coordinates has the 
drawback of messing with the well-known coordinate singular- 
ity at the poles of the spherical grid. From earlier calculations 
of the Rayleigh-Taylor mixing in core collaps e supernovae em- 
ploying spherical coordinates it is known fsee lKane et all200Ql 
and Paper I) that Rayleigh-Taylor fingers grow faster along the 
poles than in other regions of the grid, and may even evolve 
to (artificial) "axial jets" if the simulation is followed to late 
times. This eff'ect is most likely enhanced by the imposed ax- 
isymmetry, which restricts the degrees of freedom of the flow, 
and promotes its convergence along the impenetrable poles. 

If the initial conditions do not deviate too strongly from 
spherical symmetry, the axial artifacts remain rather localized, 
and do not afl'ect the global character of the solution (see the 
calculations of Paper I, and Kaneetal. 2000). However, this 
may no longer be the case if one starts from very anisotropic 
situations. To illustrate this, we show the distributions of en- 
tropy and density at ^ Is and at ^ 3000 s, respectively, for 
Models b 18a, bl8b, andb23ain Fig.[Tl Of these models the first 
one (Model bl8a) is the most anisotropic, being dominated by 
the / = 1 (bipolar) mode one second after core bounce, while 
the latter two show a dominance of the / = 2 (quadrupole) mode 
(with a smaller / = 1 contribution). 

It is apparent from Figs.[lji-f that in all of the models the 
solution at a time of ~ 3000 seconds is aff'ected significantly 
in the polar regions, showing a pair of polar "bulges" whose 
origin is in all cases a pair of polar Rayleigh-Taylor fingers. 



These expand with much higher velocity than the surrounding 
material and ultimately deform the main shock. 

Yet, while in Models bl8b and b23a the solution is still 
undisturbed at distances > 15 - 20° from the poles, the artifacts 
are very pronounced in the most anisotropic model, bl8a. In 
this case they have apparently also interacted with the physi- 
cal instability discussed in Sect. 14.21 aff'ecting a cone as large 
as ~ 45° around the south pole. We feel that the occurrence of 
a possibly spurious solution over such a large fraction of the 
computational domain makes this model unsuitable for reason- 
able analysis, and we will thus skip it from the further discus- 
sion. Hence, we will not follow the late-time evolution of this 
strongly one-sided explosion (resulting from a clearly domi- 
nant / = 1 mode) in this paper, but will constrain ourselves to 
explosions in which the / = 2 mode is dominant, and the / = 1 
mode yields only a smaller contribution. We will also try to ex- 
clude the regions of these latter models which are aff'ected by 
the "jets" from all our analyses of expansion velocities and the 
extent of mixing of diff'erent elements, by skipping 15° of the 
computational wedge closest to the north and south poles from 
the data evaluation. 

For similar reasons we will also not attempt to calculate the 
late-time evolution of rotating models in this paper. It is actu- 
ally true that in a rotating model the polar axis, i.e. the axis 
of rotation, represents a distinguished direction of the physical 
system, which is impenetrable for the flow if the specific angu- 
lar momentum of fluid elements is conserved. In this case the 
hydrodynamic solution may in fact yield the formation of jet- 
like polar outflows which are physical. In 2D (i.e. axisymmet- 
ric) calculations, employing standard spherical coordinates, it 
is, however, impossible to study the formation of such features 
reliably, since the enforced conservation of the z-component 
of the specific angular momentum leads to numerical artifacts 
near the poles which will always contaminate the numerical 
solution. 

Unfortunately, a reliable modelling of late-time supernova 
hydrodynamic instabilities (beyond times of ~ 10 s) for explo- 
sions with extreme initial anisotropics and/or rotation appears 
to require high-resolution 3D simulations of the entire sphere, 
either in Cartesian coordinates, or with a compo site spherical 
mesh, as e.g. the "cubed sphere grid" described in lRonchi et al.l 
(Il996l) . Both approaches are expected to result in computa- 
tional tasks that are between two and three orders of magnitude 
more expensive than the present calculations. 

4. Results 

4.1. Energetics, timescales, and morphology 

Table |2| provides an overview of all models that we have 
computed for the 15 blue supergiant progenitor of 
Woosle v et^ (Il988h to date. Model bl8a (which we wiU not 
discuss further due to the reasons mentioned in Sect. 13.21) and 
Model bl8b are "medium core luminosity" models with an ex- 
plosion energy of 10^^ erg. They diff'er only in the random seed 
perturbation that we added to the velocity field of Bruenn's 
post-bounce model to trigger the onset of neutrino-driven con- 
vection. They diff'er hardly in their global parameters like the 
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Fig. 2. Snapshots of the density distribution of Model b23a for several times. From top left to bottom right a) ^ = 4 s, b) ^ = 10 s, 
c) r = 20 s, d) ^ = 50 s, e) ^ = 99 s, and f) t = 299 s. Note the change of the radial scale. The supernova shock is the outermost 
discontinuity which becomes progressively spherical with time. Note also the growth and outward propagation of the Rayleigh- 
Taylor mushrooms from the former Si/0 and 0/He interfaces of the star, and the onset of the Richtmyer-Meshkov instability at 
the He/H interface (the deformed discontinuity just behind the supernova shock in the plots for t = 99s and t = 299 s). 
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Fig. 3. Same as Fig.|2l but for selected late instants in time. From left to right a) t ^ 1500 s, b) t ^ 10000 s, and c) t ^ 20000 s. 
The outer circular boundary in the last two plots is the outer boundary of the computational domain, which the supernova shock 
has already left at those times. The shock is only visible in the plot for t ^ 1500 s and has become almost completely spherical 
at this time, while it is decelerating in the hydrogen envelope. Below it, the Richtmyer-Meshkov instability at the He/H interface 
has grown, while in still deeper layers a strong reverse shock decelerates the matter of the He core. Note (by comparing Fig.|3t 
with Fig. |2f) that a few of the Rayleigh-Taylor clumps have actually reached the He/H interface before this reverse shock has 
managed to form, and that thereby they have escaped an interaction with this shock. Note also the huge hydrogen pockets that 
are created by the Richtmyer-Meshkov instability. 



Table 2. Models for the^Woos lev et al.l 11988*) blue supergiant. 
Given are the explosion timescale, ^exp, and the entire simulated 
time, ^sim, as well as the explosion energy at a time of 3000 
seconds, E^^^^^^^, and the baryonic mass of the neutron star at 
t ^ I s (without fallback). The explosion timescale is defined 
by the time after core bounce when the integral of the total 
energy over all grid cells with positive values of that quantity 
and positive radial velocities, is larger than 10"^^ erg. 



Model 


^exp 


^^=30008 




^sim 




[s] 


[10^^ erg] 


[Mo] 


[s] 


bl8a (disregarded) 


0.190 


1.0 


1.3 


10^ 


bl8b 


0.185 


1.0 


1.3 


2x 10^ 


b23a 


0.138 


2.0 


1.2 


2x 10^ 


T310a (Paper I) 


0.062 


1.7 


1.1 


2x 10^ 



explosion energy, but show large diff'erences in the morphol- 
ogy, indicating the extremely nonlinear character of the growth 
of non-radial instabilities in the SN core. Model b23a, on the 
other hand, is a "high core luminosity", high-energy explosion 



of 2 X 10^^ erg. We recall once more that in all of these sim- 
ulations the core luminosities were held constant for the first 
second after core bounce (see Sect. 13. 11) . 

In contrast. Model T310a (which has been described in de- 
tail in Paper I) employed core luminosities that declined ex- 
ponentially with time. We will include it in the following dis- 
cussion in order to demonstrate that our new calculations show 
dramatic quantitative and qualitative diff'erences compared to 
this older simulation. We strongly encourage the reader to com- 
pare the density and velocity plots of Model T3 10a that we pre- 
sented in Paper I, with the plots for Model b23a (Figs.|2l[3land 
and Model bl8b (Figs. 0] and that we show in this paper. 
The diff'erent morphology of the neutrino heated ejecta and the 
shock within the first few seconds after core bounce is striking. 
As already stated in the introduction, this is basically a conse- 
quence of diff'erent explosion timescales and the thus diff'erent 
growth of non-radial instabilities during the shock revival phase 
by neutrino heating, and the onset of the explosion. 

With a timescale ^exp of only 62 ms. Model T3 10a exploded 
fast - nearly as fast as the typical turnover time of a convec- 
tive eddy. This neither allowed convection nor the advective- 
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Fig. 4. Snapshots of the density distribution of Model bl8b for several times. From top left to bottom right sl) t = 4s,h) t = 100 s, 
c) ^ = 299 s, d) ^ = 1512 s, e) ^ 10 000 s, and f) t ^ 20 000 s. Note the change of the radial scale. The supernova shock is visible 
in plots a)-d) only. The outer circular boundary in plot f) is again the outer boundary of the computational domain. Note also the 
more pronounced hemispheric asymmetry of the ejecta at late times as compared to Model b23a shown in Fig.fS] 
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acoustic shock instability to grow appreciably. Only small con- 
vective blobs were formed in the neutrino-heated layer out 
of the imposed random perturbations, and the shock wave re- 
mained undistorted (i.e. spherical). Moreover the flow pattern 
was dominated by a high-order mode, with about a dozen ed- 
dies. 

In contrast. Models b23a and bl8b exploded more slowly 
by a factor of between two and three (Table |2|). Hydrodynamic 
instabilities in the post-sh ock flow thus had more time to grow. 
As in the simulations of IScheck etakl (I2QQ4 l2QQ6h the ini- 
tially small blobs in the neutrino-heated layers merged, and 
ultimately formed a low-order, 1 = 2 dominated mode with 
a smaller 1=1 mode contribution, exhibiting only two large 
buoyant bubbles which are separated by a single accretion fun- 
nel. In the process, they deformed the shock wave to a prolate, 
almost peanut-like shape (see Figs.[lJ),[l]:,|2t, and|4t). 

The longer explosion time scales also led to larger neutron 
star masses in Models b23a and bl8b around one second af- 
ter core bounce than in Model T310a (Table |2|). While these 
masses are still rather low, the comparison shows that the exact 
value is fairly sensitive to the treatment of the neutrino physics 
and the assumed boundary conditions (i.e. the neutrino prop- 
erties and the contraction behaviour of the neutron star core). 
When the explosion sets in later, the neutron star mass is higher. 
Moreover it also depends on the structure of the progenitor star 
(see Scheck et al. 2 006 for the corresponding significant varia- 
tions in case of diff'erent 15 M© progenitors), and can increase 
due to later fallback. We have not traced this fallback in the 
present work, though. 

4.2. Shock expansion and Richtmyer-Meshkov 
instability 

The early global anisotropy of the shock in the low-mode mod- 
els is decisive for the further evolution. It leads to qualita- 
tively new efl'ects in the late-time hydrodynamics compared to 
Model T310a, in which the shock rema ined spherical f r om th e 
outset, or the simulation presented in Kifo nidis et al.l (|2P00^, 
where the shock had been deformed by rising bubbles during 
the neutrino-heating phase, but had become spherical by the 
time it emerged from the He core (see ' Chevalier & SokeJl989; 
for a discussion of the underlying mechanism). Note that in 
our previous simulations such deformations of the shock due to 
buoyant neutrino-heated bubbles were in all cases local ones, 
with deviations of the local to the mean shock radius of at most 
20%. In contrast. Model b23a has a prolate shock with an ini- 
tial major to minor axis ratio of ~ 1.5. A deformation of this 
magnitude not only makes the definition of a mean shock ra- 
dius questionable. It is so substantial that it apparently cannot 
be smoothed out completely by the time the blast arrives at the 
He/H interface of the Woosley et al. progenitor. 

In Model b23a the original peanut-like shock gives way to 
a surface with an equatorial bulge, that forms at the "peanuts' 
waist" a few seconds after core bounce, and progressively at- 
tempts to bring the shock into the spherical shape (Fig.|2|l. Yet, 
it requires some time to succeed in this, and when the blast has 
crossed the He/H interface around 100 s after core bounce its 
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Fig. 5. Integrated vorticity Q according to Eq. Q as a function 
of time during the late evolutionary phases of Model T310a 
(dotted line). Model b23a (solid line), and Model blSb (dashed 
line). Note that Q is essentially zero for the high-mode model 
T310a, while it is of the order lO^^cm^s"^ for the low-mode 
models b23a and bl8b, which experience a strong Richtmyer- 
Meshkov instability at the He/H interface. 



surface still shows two kinks. Kinks in the shock surface are 
usually found above prominent downflows, and our "two-kink 
feature" originated from the merging of two such downflows. 

The kinks in Model b23a are located at similar latitudes in 
the northern and southern hemispheres (see Fig.|2t). In these 
spots the shock hits the composition interface obliquely. As a 
result, it deposits a significant amount of vorticity into the inter- 
face layer, which triggers strong Richtmyer-Meshkov mixing 
(|Richtmyer 1960; Meshkov 1969), much as in the simp l e, pla- 
nar geometry setup considered by iHawlev & Zabuskvl (Il989l) 
(in fact, our Fig. [5^ shows all features that are visible in the 
plots of their flows, in particular of the "fast- slow" setup that 
they studied). Two prominent vortices^ form, which transport 
hydrogen-rich gas in the form of two large "pockets" into the 
helium and metal core. At the same time material from these 
core layers is dragged outward into the hydrogen envelope, 
while a vortex sheet due to the Kelvin-Helmholtz instability 
forms along a significant part of the interface (Figs.[3t-c). 

Note that the He/H interface is actually susceptible to both, 
the Richtmyer-Meshkov, and the Rayleigh-Taylor instability, 
which is due to pressure and density gradients of opposite 
sign in this region (see Paper I). At this interface, however, 
Richtmyer-Meshkov instabilities grow faster and appear earlier 
than RT mushrooms^ (compare Figs.|3j) and [31:). On the long 
run the mixing in fact is a result of both instabilities interact- 
ing with each other. Yet, we will still refer to the mixing at the 

^ In the context of axisymmetric simulations all flow features are 
actually toroidal structures, and the "vortices" are actually "vortex 
rings". 

^ The Richtmyer-Meshkov instability is an essentially impulsive in- 
stability triggered upon shock passage. The RT instability is taking 
place on a much longer timescale due to the smooth large-scale gradi- 
ents of pressure and density at the He/H interface. 
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He/H interface (a bit loosely) as "Richtmyer-Meshkov instabil- 
ity" in order to make clear that it is the deposition of vorticity 
by the aspherical shock which is crucial. To demonstrate this, 
we consider 

n = I (V X w) dV, (3) 



the integral of the vorticity over the computational domain. In 
Eq. ^ u is the flow velocity, V is the computational volume, 
and dV is the volume element. For 2D axisymmetric calcula- 
tions in spherical coordinates there is only one component of 
V X M that does not vanish, namely (V x m)^, and therefore Q 
reduces to the scalar 



r^sin?^drd?^, (4) 



where u and v are the velocity components in r and direction, 
and Ri and Rq are the inner and outer (open) radial boundaries, 
respectively. Note that Q is a suitable measure for the strength 
of vortices residing in the domain. 

Figure 13 shows the temporal evolution of Q for Models 
b23a, bl8b, and T310a, and times between 10"^ and 2 x 10^ s. 
When the shock has already left the grid (which is the case at 
these times), no contamination of Q due to the axial artifacts 
visible in Fig. [I] is present any more, and the only significant 
contribution stems from the vorticity deposited early on by the 
shock at the He/H interface. This contribution is of the order of 
a few 10^^ cm^s"^ for Models b23a and bl8b, while it is essen- 
tially zero for Model T310a with its spherical shock. 

We have already shown in Paper I that Model T3 10a exhib- 
ited only a very weak Rayleigh-Taylor instability at the He/H 
interface, which was triggered by perturbations of this interface 
due to acoustic noise (i.e. sound waves). Although this led to 
the growth of a large number of small Rayleigh-Taylor mush- 
rooms. Model T310a clearly /a/Z^J to mix significant amounts 
of hydrogen into the helium and metal core, i.e. below a mass 
coordinate of ~ 4Mq (compare Fig.|6t). The "high- vorticity" 
models, b23a and bl8b, on the other hand, succeed in mixing 
hydrogen down to a mass coordinate of ~ 1.3 M© (Figs. [St and 
c), which is close to the neutron star mass. This clearly demon- 
strates that in addition to the short- wavelength Rayleigh-Taylor 
instability (which can be easily identified in the late evolution- 
ary stages of Models b23a and b 18b by the typical mushroom- 
shaped structures, see Figs. 130] and Q, a much more efficient 
vortical mixing mechanism is present in these models, which 
operates on much larger scales. Here we identify this vortical 
mechanism with the Rich tmyer-Meshkov instability (see also 
iHawlev & Zabuskvl 19891 and the references therein). 

A consequence of the nature of this instability is that it 
leads to a dependence of the final (large-scale) spatial distribu- 
tion of chemical elements on the initial shape of the shock, i.e. 
on the spectrum of the superposed unstable modes in the early 
post-shock flow pattern. If, initially, two equally large buoyant 
high-entropy bubbles are present behind the shock, which are 
separated by a single prominent accretion funnel that is located 
close to the equatorial plane, i.e. if a dominant / = 2 mode is 
present (see Figs, [l]: and |2^), then the later distortion of the 
shock surface still reflects this initial asymmetry of quadrupo- 
lar character. In this case two nearly equally large pockets of 
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Fig. 6. Extent of the chemical mixing in the simulations at a 
time of 10 000 s. From top to bottom a) Model T3 10a, b) Model 
b23a, and c) Model bl8b. In all cases the angle-averaged mass 
fractions of H (bold soHd line), "^He (bold dotted line), ^^O (thin 
dashed line), ^^Si (thin dotted line), ^^Ti (bold dashed-dotted 
line), and the iron group (thin solid line), are shown as a func- 
tion of the enclosed mass. Note the deep inward mixing of hy- 
drogen and the strong outward mixing of the metals in Models 
b23a and bl8b, and the lack thereof in Model T310a. 



hydrogen finally develop, which, together with the distribution 
of the outward dragged material from the helium and metal 
core, show almost equatorial symmetry (Fig.Et)- If on the other 
hand, the buoyant bubbles are of diff'erent size and the accre- 



12 



K. Kifonidis et al.: Non-spherical core collapse supernovae 



tion funnel is bent and inclined towards one of the poles, i.e. if 
there is some admixture of the 1=1 mode in the initial flow 
pattern, as in Model bl8b (see Figs.[TJ) and0t). the H pockets 
tend to grow to difl'erent sizes. The final hydrogen distribution 
then develops a pronounced hemispheric asymmetry, as do also 
the spatial distributions of helium and the metals (Fig.Cb). 

Both the deep inward mixing of hydrogen, and the natu- 
ral development of large-scale, hemispheric asymmetries in the 
spatial distribution of diff'erent elements, have important con- 
sequences for the interpretation of observational data of su- 
pernova explosions and supernova remnants. The former is re- 
quired to obtain good fits in light curve modelling, while the 
latter is probably the cause of the asymmetric iron (and nickel) 
lines of SN 1987 A. We will return to these issues in Sect. 13 

It is quite astonishing that (unlike the RT instability) the 
importance of the Richtmyer-Meshkov instability for the obser- 
vational appearance of core collapse supernova explosions has 
not been widely recognized in previous modelling (but see e.g. 
iKan e et al .'2000 for an exception). This may be, at least in part, 
due to the fact that in previous studies of aspherical shock prop- 
agation through exploding stars this instability was not noticed. 
Judging from the published infor mation, it is fo r instance not 
visible in the calculations of Hunger ford et al.l (2003, 2005). 
These authors started their 3D simulations at a time of 100 s 
after core bounce from a ID calculation of the earlier phases 
of the explosion. Thus the shock had already crossed the He/H 
interface at the time they added asymmetries to the velocity 
field of their spherically symmetric explosion models. This ob- 
viously gave the instability no chance to develop. 

In the works of |Na^ataki| §000) and |Yamada&Sat3 
(Il99lh . on the other hand, who initiated their 2D simulations 
from aspherical shocks at much earlier times, it was proba- 
bly the coarse angular resolution of only 100 zones which pre- 
vented them from discovering the instability (in the latter paper 
a diff'usive advection scheme with only first-order accuracy in 
space was employed, too). We think that this demonstrates the 
importance of two points which we have already addressed in 
Paper I, namely that modelling the mixing in core-collapse SNe 
requires one to follow the evolution from the earliest phases of 
the explosion in more than one dimension (in order to avoid 
biasing the late-time evolution). Furthermore, high resolution 
and an excellent advection scheme are essential for resolving 
and properly describing the relevant fluid instabilities. 

4.3. Metal clumps and Rayleigh-Taylor instability 

The unstable modes that become dominant within the first sec- 
ond of the explosion, do not only determine the morphology 
of the shock, and thereby the sequence of events at the He/H 
interface. They also provide the perturbations for the Rayleigh- 
Taylor unstable Si/0 and 0/He interfaces of the progenitor. We 
have noticed in Paper I that the wavelength of these pertur- 
bations determines the number of Rayleigh-Taylor fingers that 
grow at these interfaces, and that it characterizes the flow pat- 
te rn until late ti mes. This is also supported by the recent work 
of 'Miles ('2004'), who has developed a Rayleigh-Taylor bubble 
merger model for compressible flows. He concluded that un- 



less the shock Mach number and the mode number of the initial 
perturbation are very high, the system will not loose memory 
of the original perturbation spectrum. 

This is indeed the case in all models that we have computed 
to date. Model T3 10a developed about a dozen Rayleigh-Taylor 
fingers out of a comparable number of flow features that formed 
during the phase of neutrino-driven convection near the Si/0 
and 0/He interfaces (Paper I). The low-mode models b23a and 
blSb, on the other hand, lead to the growth of only about three 
to four, rather large, Rayleigh-Taylor mushrooms (Figs. |2| and 
|4]|. In what follows we will keep referring to the case of Model 
T310a as a high-order, and to this of Models b23a and bl8b as 
a low-o rder mode, in order to distinguish them clearly. In the 
sense of lMilesI (l2004l) all of our models bear "low-mode" per- 
turbations, because they all satisfy his condition for retainment 
of the initial mode number. 

Besides having a difl'erent number and size of Rayleigh- 
Taylor "clumps", both cases also difl'er strongly in the velocity 
distributions of nucleo synthetic products. The low-mode mod- 
els show much broader distributions in velocity space for iron- 
group nuclei and for the isotopes ^^O, ^^Si and ^^Ti a few sec- 
onds after bounce than Model T310a (compare the plots for 
t = 4sm Figs. [8land|9l with Fig. 16 of Paper I). 

Even more importantly, the maximum velocities of these 
nuclei are significantly larger than in the high-mode model 
T310a. In Model b23a we find maximum metal velocities 4 
seconds after bounce that are about 40% larger than in Model 
T310a. But even in the low-energy, low-mode model blSb, 
whose explosion energy is nearly half that of Model T3 10a (1.0 
vs. 1.7x10^^ erg, respectively), the corresponding velocities are 
by about 20% larger. 

The spatial distribution of kinetic energy in the low and 
high-mode cases is fundamentally diff'erent. In Model T310a 
most of the kinetic energy was contained in an essentially 
spherical shell behind the shock (similar to a ID simula- 
tion), and there was only little kinetic energy in the two- 
dimensionally perturbed region containing the bubbles. The 
weak neutrino-driven convection in Model T310a was not able 
to boost the velocities inside those bubbles relative to the mean 
expansion velocity of the background flow. In the low-mode 
models this is difl'erent. The high-entropy bubbles and the 
higher-density pockets enclosed by them contain most of the 
kinetic energy of the explosion. 

The higher metal velocities during the first few seconds 
after core bounce lead to significant diff'erences in the late- 
time evolution of the low-mode models as compared to Model 
T310a. The metal clumps now propagate faster through the He 
core of the star, keeping closer contact to the outward sweep- 
ing supernova shock. In both. Model b23a and Model bl8b, the 
fastest clumps have almost reached the He/H interface around 
300 s after core bounce (see Figs. |2f and |4l:). At this time the 
strong reverse shock that will develop below this interface due 
to the deceleration of the main shock in the hydrogen envelope 
(compare Fig. [3^ and Paper I), has not yet formed. Hence, an 
interaction of the fastest clumps with this reverse shock, which 
leads to the dramatic slow-down of the metals beyond a time 
of 1500 s in Model T310a (Paper I), does not happen in the 
low-mode models ! 
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Fig. 7. Entropy distribution (in units of ^B/nucleon) of the low-mode models at late times. From top to bottom: a) Model b23a at 
^ 10 000 s and b) Model bl8b at ^ 20 000 s. The color coding has been chosen such that material with different composition 
can be clearly distinguished. Hydrogen-rich gas is shown in orange or red, while material of the helium core is coded in dark 
green. Matter from the metal core, which includes nuclei that have been synthesized during the explosion as well as during the 
hydrostatic evolution of the star, is shown in light and dark blue. 



Instead, the fastest clumps encounter a layer in which the 
density gradient, though steepening, is still relatively flat. Their 
motion relative to the background flow thus never becomes su- 
personic, and deceleration enhanced by supersonic drag is ab- 
sent. Hence the maximum velocities of iron-group nuclei and 
other elements decrease only slightly at late times, e.g. from 
3500 km/s at t = 1500 s to 3300 km/s at r = 10 000 s in Model 
b23a (Fig.[8l). These velocities remain unchanged until the time 
when the clumps start to leave the outer boundary of our grid, 
which happens at ^ 12 000 s in Model b23a. 

Strong reduction of a late-time deceleration may also be 
favoured by the fact that the metal clumps become engulfed 
in the Richtmyer-Meshkov instability at the He/H interface. 
The vortical energy deposited early on in these layers by the 
shock, represents a reservoir of kinetic energy that can pos- 
sibly be tapped during the late-time hydrodynamic evolution, 
and would have been unavailable to the clumps if the shock 



had been spherically symmetric when it crossed the interface 
(as it was the case in Model T310a). 



Summarizing these findings, we can state that by following 
the late-time hydrodynamic evolution of neutrino-driven explo- 
sions which were dominated by early, low-mode hydrodynamic 
instabilities, we have for the first time obtained a consistent (in 
the sense defined in Sect. [B solution to the nagging "nickel 
discrepancy" problem (Hera nt & Benzlfl992l) . i.e. we have ob- 
tained maximum velocities for the bulk of the metal-rich ma- 
terial (in Model b23a) which are in good agreement with those 
of SN 1987 A. This result calls for a more careful inspection of 
the simulations in order to judge whether they also fulfill other 
observational constraints obtained from this supernova. 
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Fig. 8. Logarithm of the fractional mass of different elements as a function of the radial velocity Vr for Model b23a at various 
epochs. Note the different scales for Vr in the left and right panels. Note also that the fractional mass was computed with respect 
to the mass of the element, M/, that is contained on the current extent of the computational grid. For the lighter nuclei H and 
^^O, this mass is not necessarily identical to the total element mass, because initially our adaptive mesh does not yet cover the 
outermost stellar layers, and because in the late phases the shocked gas is allowed to stream off the grid. The latter fact is the 
reason why the hydrogen velocity distribution at ^ = 10000 s does not extend beyond ~ 4000 km/s (because faster expanding 
hydrogen has left the grid already). 



5. Comparison with SN 1987 A 

The most satisfactory approach for verifying whether our simu- 
lations are consistent with the observational data of SN 1987 A 
would, of course, be the computation of light curves, spectra 



and the degree of polarization from the 2D hydrodynamic mod- 
els. This is beyond the scope of the present paper and will be 
attempted in future work. 

Here we will confine ourselves to a semi-quantitative com- 
parison with a check list that we assembled from the results of 
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Fig. 9. Same as Fig. El but for Model bl8b. 



previous observational and (one-dimensional) modelling work. 
This list, which any successful model of SN 1987 A needs to 
pass, consists of the following items: 

1 . The high velocities of iron-group nuclei, or, given that ul- 
timately homologous expansion establishes, the outward 
mixing of these nuclei in mass. We will focus here on the 
mixing in velocity space since this could be measured di- 
rectly by observations of the Fe and Ni infrared lines. The 
mixing in the mass coordinate can be deduced only indi- 
rectly, on the basis of theoretical models, and may thus be 
less reliable. 



iHaas et al.l (I199Q|) and lColgan et al.l (Il994 have carried out 
observations of the [Fe II] 18 /im, [Fe II] 26 yum, and [Ni I] 
6.6 yum lines at 410 and 640 days after the explosion. All of 
these lines show broad wings. In addition, there is an iso- 
lated high- velocity emission feature in the red wing of both 
the iron lines at day 410, and the nickel line at day 640, 
which indicates that a single clump with ~ 3% of the total 
iron mass mo ves away from the observer with a vel ocity 
of 3900 km/s ( Haas et al. 1990). Jennings et al.' ("1993) give 
an expansion velocity of 3150 km/s for this high- velocity 
feature, based on independent observations of the [Fe II] 
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18 um line at day 6 13. The [Ni I] 6.6 fim measurements of 
IColganetal.1 (11994 at day 640 exhibited probably the high- 
est signal-to-noise ratio but were affected by absorption due 
to dust, which had for med in the ejecta between days 410 
and 640. Colgan et all (1 1994 thus regard the [Fe II] lines 
at day 410 as providing the best evidence for the amount 
of iron-group elements at high velocities. From the width 
of these lines Haas et al. (1990) concluded that > 4% of 
the total iron (excluding the high- velocity feature) has an 
expansion velocity > 3000 km/s. 

Early nickel mixing with significantly hi gher velocities 
(above 5000 km/s) has been postulated bv ^ Mitchell et al.l 
(l200lh in order to get sufficient excitation for the hydrogen 
B aimer lines in their steady state spectral models for day 
4. IUtrobin & Chu gar(2005). however, have recently shown 
that such high (bulk) nickel velocities are not required to 
explain the observations during the photospheric phase, if 
one takes into account the effects of time-dependent ion- 
ization in modelling the spectra. With their time-dependent 
ionization code excellent fits to the Ha profiles are ob- 
tained if the ^^Ni is mixed out to only ~ 3000 km/s. Further 
work which provides compelling evidence of such more 
moderate velocities for the bulk of the ^^Ni (and, for in- 
stance, includes modelling of the X-ray and r-r ay emis- 
sion of SN 1987 A) is summarized in the review o f lMcCravl 
(Il993h. So me more recent papers are the light curve study 
of Utrobin (2004), and the work of F assia & M eikle ( 1999). 
The latter authors have applied spectral synthesis models to 
the He 1 10830 A line after day 10 and found that about 3% 
of the total ^^Ni mass must be mixed to above 3000 km/s, 
but that the ^^Ni concentration must be negligible above 
3900 km/s. They also demonstrated that models with this 
amount of mixing are able to reproduce the y-ray light 
curves and the late time infrared line velocities. 
The marked asymmetry of the [Fe 111 and [Ni I] in- 
frared fines (Haasetal. "l990"; Snvro milio et al.l [1990; 
IJenningsetak 1993: Cokan et al. 1994) . Besides the (lo- 
calized) "high- velocity-feature" that we already mentioned 
above, these lines also show a global red/blue asymmetry. 
Its origin was pro bably anisotropic expansion of the eje cta 
on a global scale JHaas et alll99Ql:ljennings et alJI 19931) . 
The total ^^Ni mass of ~ 0.07 Mo, which is required to 
power the li^ht curve (^ ,W ooslev 1988: Suntzeff' & Bouchet 



3. 



BlBouchetet alii 991 

4. The very low minimum hydrogen vel ocities. The mere fact 
that the Ha line profile observed bv IPhillips et al. (,1990) 
on day 498 is not fiat- topped implies that there is no large 
hydrogen cavity in the slowly expanding central layers of 
the ejecta (V. Utrobin, private comm unication). From the 
narrow core of t he Ha line profiles of IPhillips et al . ( 1990) 
around day 800.'Koz ma & Fransson ! (^1998') inferred a min- 
imum hydrogen velocity of < 700 km/s. Mixing of hy- 
drogen into the SN core, combined with the outward mix- 
ing of ^^Ni, is also necessary to fit the broad peak of 
the light curve (e.^ . Woosley 1988; Shigevama & Nomotd 
ll990l:lBlinnikovlll999l) . In the modelsof iBlinnikovl (1 1999) 
hydrogen is mixed down to a mass coordinate of ~ 1.3Mo. 



5. The prolate deformation of the (inner) ejecta with an axis 
ratio of roughly 2:1 that has been found by high resolu- 
tion imaging with the HST, and which has been observed 
indirectly by earli er spectropolarimetry measurements (see 
I Wang et al]l2002L and the references therein). It is partic- 
ularly noteworthy that the degree of polarization was only 
0.1% two days after the explosion and increased to 1.1% 
at about day 200, indicating that the scattering surface be- 
came increasingly asymmetric as the o bservations probed 
the deeper layers of the ejecta dWang et al...2002) . 

6. The constraints on the ratio of the explosion en- 
ergy to ejected mass, Eq^^^/M, that have been deduced 
from modelling of the light curve and the hydrogen 
lines. Utrobin & Chugail |2005| ) give a value of 0.83 x 
Iq50 gj,g M~^. Blinnikovl (^1999) finds, from a multigroup 
radiation hydrodynamic study, Eq^^/M - (0.75±0.17)x 
10^^ erg (based on ^exp = (1.10 ± 0.25) x 10^^ erg, 
and an assumed ejected mass of 14.7 M^), while from the 
models listed in Table 2 of^Shi gevama & Nomoto ( 1990) an 
average value of about 0.76 x 10^^ erg can be inferred. 



We have already noted that in the high-energy model, b23a, 
the (final) maximum iron-group velocities are about 3300 km/s. 
Furthermore, a few per cent of the total iron-group mass in 
this model is expanding with velocities > 3 000 km/s (compar e 
Fig. [HI), as required by the observations of iHaas et alJ (ll99Qh . 
Our simulations are thus able to reproduce the maximum ve- 
locities of the bulk of the iron-group mass in S N 1987 A. The 
explanation of single, even faster clumps ( Utrobin etalJll995l) 
might require different initial conditions, or higher resolution, 
or will have to await future 3D simulations. Note also that in 
our low-energy model, bl8b, the highest iron-group velocities 
are ~ 2200 km/s, which is too low and indicates that the explo- 
sion energy was higher than 10^^ erg (see also the comments 
with respect to item six below). 

To be consistent with the second item of the above list, the 
simulations are required to exhibit variations of the iron veloc- 
ities, and/or the iron concentration with respect to the latitudi- 
nal angle on a large scale, preferably in the form of a hemi- 
spheric asymmetry. In this case one would measure higher ve- 
locities, and/or a larger line-flux from one hemisphere of the 
ejecta than from the other. Our models indeed develop such 
hemispheric asphericities, if a sufl&ciently large contribution 
of the / = 1 mode to the shock morphology is present after 
the neutrino-heating phase. This, on the other hand, depends 
on the imposed initial perturbations, which afl'ect the forma- 
tion and merging of the neutrino-heated bubbles during the first 
~ 300 ms of the explosion. S ince this merging is a chaotic and 
highly nonlinear process (^cf. lScheck et alJl200^ . the develop- 
ment of a hemispheric asymmetry is in the end governed by 
stochastics. Model bl8b shows such a pronounced asymmetry, 
while in Model b23a (which fits most of the other check- list 
items better than Model bl8b) such an asymmetry is there, too, 
but weaker (Fig.Q). Yet, the outcome might have well been re- 
versed if different initial perturbations were used. We refrain 
from demonstrating this here, as this would require us to calcu- 
late a large sample of models. 
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A reproduction of item three hinges very sensitively on the 
electron fraction, Yq, of the neutrino-heated ejecta. The local 
Ye is in turn a result of the competition of the charged-current 
reactions 

Vq -\- n ^ e~ -\- p , (5) 
Vq -\- p ^ e'^ -\- n . (6) 

Electron capture on protons and Vq absorption by protons de- 
crease Fe and make the matter more neutron rich. On the other 
hand, the inverse processes of positron capture on neutrons and 
Ve absorption by neutrons drive the matter to the proton-rich 
side. The rates for the Ve and Ve captures depend sensitively on 
the local Ve and Ve spectra, and hence they require an energy- 
dependent transport treatment for their accurate calculation. 
With a gray transport scheme like ours this can be done only 
approximately, resulting in an uncertainty in Yq during nuclear 
freezeout of the order of several per cent. Since a variation of 
only a few per cent in Yq can result in a more than 100% change 
in the ratio of the ^^Ni yield relative to the yields of other (e.g. 
more neutron rich) iron-group nuclei, it is not possible from the 
present calculations to determine reliable isotopic yield ratios 
within the iron-group. We can, however, give an upper limit 
for the ^^Ni yield, obtained from the Fe -insensitive total yield 
of iron-group nuclei. This upper limit amounts to 0.15 M© and 
0.11 Mq in Models b23a and bl8b, respectively, and is com- 
patible with the fact that 0.07 M© of ^^Ni were produced in SN 
1987 A. 

The consistency of the simulations with item four is easily 
verified from Figs.[8land |9l which show that the lowest hydro- 
gen velocities in both Models b23a and blSb are ~ 500km/s. 
That the models are also consistent with the extent of hydro- 
gen mixing assumed in lightcurve modelling is best checked by 
comparing our Fig. |6l with Fig. 2 of Blinnikov ( 1999). The hy- 
drogen mass fraction profile used bv lBlinnikovl to obtain a good 
fit to the observed lightcurve is almost indistinguishable from 
the ones that we obtain for Models b23a and bl8b. The sharp 
drop interior to ~ 1.5 M© and even the plateau at logX = -1 
for mass coordinates 1.5 M© < Mr < 3.5 M© are reproduced 
excellently. 

Consistency with item five can be inspected qualitatively 
from Figs. OH andQ which demonstrate that at late times the 
inner ejecta (i.e. the helium and metal-rich material) in our sim- 
ulations exhibit a global anisotropy. The major to minor axis 
ratio of the isodensity contour corresponding to the He/H inter- 
face is roughly 1.6. The decrease of the asymmetry of the ejecta 
with radius inferred from the spectropolarimetry data might be 
explained from the fact that the shock in our models (baring the 
numerical artifacts near the poles) becomes spherically sym- 
metric once it has propagated through the inner hydrogen en- 
velope (see the plots for 1500 s in Figs.[3land |4|). This is an im- 
portant feature of our models which clearly distinguishes them 
from earl ier simulations, e.g. fro m the asymmetric models pre- 
sented in Yamada & Sato ( 1991). A more detailed comparison 
with the spectropolarimetry data - beyond the qualitative level 
that we attempt here - is difficult because it requires a calcula- 
tion of the (time-dependent) location of the photosphere in our 
models, which we did not perform. 



This leaves us with item six. For our 15 M© progenitor and 
Model b23a, with its explosion energy of 2.0 x 10^^ erg, we 
obtain an energy-to-mass ratio of 1.4 x 10^^ erg (after ac- 
counting for a neutron star of 1 .2 M© and ignoring possible fall- 
back). This is a bit on the high side. For the low-energy model 
blSb we get 0.7 x 10^^ erg , which is in much better agree- 
ment with the lightcurve modelling, but for this model the max- 
imum iron-group velocities are only ~ 2200 km/s, which is in 
conflict with item one. Still, we think that with a more careful 
choice of the progenitor model and of the explosion energy, it 
will be possible to obtain both, agreement with the observed 
maximum iron-group velocities, and a closer reproduction of 
the £'exp/M ratio because 

- the mass of the SN 1987 A progenitor may have been in the 
18 - 20 M© range (Wooslev 1988), 

- a lowering of the explosion energy from 2.0 x 10^^ erg to 
1.5x10^^ erg will lower our maximum metal velocities only 
by a factor of ~ 1.15 (because these velocities scale approx- 
imately with the square root of the explosion energy), 

- in three dimensions it is expected that, for the same explo- 
sion energy, larger maximum metal velocities will be ob- 
tained (due to the small er drag experience d by truly three- 
dimensional clumps, seelKane et allbOOOl). 

In summary, the low-mode models, and in particular the 
high-energy case b23a, are in at least qualitative agreement 
with the items listed above. In many cases even an astonish- 
ingly good quantitative agreement has been obtained. This is 
remarkable, especially if one considers the fact that there is 
only very little room for "fine-tuning" of parameters in the cal- 
culations that we have presented here. The freedom we have 
to adjust the model is limited to the choice of the initial ran- 
dom seed perturbation (which influences the relative contribu- 
tion of different low-order unstable modes in the post-shock 
ffow during the shock revival phase), and the parametrization 
of the inner boundary condition which determines the (elec- 
tron ffavour) neutrino ffux from the neutron-star core, L^°^^, and 
the reby the explo sion time scale and energy of the simulation 
rcf. lScheck eta TI EoO^ . Note also that the core neutrino ffux 
is only a part of the total neutrino ffux in our calculations. A 
signiffcant part of the total neutrino luminosity is contributed 
by the neutrino emission of mass accreted onto the neutron 
star, which is taken into account by our transport scheme (see 
IScheck etalJl2006l) . The consistency of our models with the 
above list therefore makes us conffdent that the long-standing 
mystery concerning the origin of the major observational prop- 
erties of SN 1987 A may have found its solution. 

6. Conclusions 

The main conclusion of the present study is that if a low- 
mode hydrodynamic instability (more precisely, a quadrupole, 
1 = 2, mode-dominated instability with a smaller contribution 
from a dipole, 1=1, mode) establishes during the ffrst second 
of a supernova in a stan dard 15 M© blue supergiant star (see 
Scheck et al. 2004, 2006 for the requirements of this), then the 
ensuing explosion possesses properties which are in remark- 
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able quantitative agreement with the observational data of SN 
1987 A. 

Among other features, explosions of this type (with an en- 
ergy of ~ 2 X 10^^ erg for the progenitor model we used) exhibit 
final iron-group velocities of up to 3300 km/s, strong mixing at 
the He/H composition interface, with hydrogen being mixed 
downward in velocity space to only 500 km/s, and a final pro- 
late anisotropy of the ejecta with a major to minor axis ratio of 
about 1.6. 

The success of low-mode explosions to reproduce these ob- 
servations is based on two eff'ects: The long- wavelength pattern 
with which the Rayleigh-Taylor unstable Si/0 and 0/He inter- 
faces of the presupernova are perturbed, and the initial global 
deformation of the shock, which induces a strong Richtmyer- 
Meshkov instability at the He/H interface^ 

Compared to our previous models jKifonidis et aLlEoQ^I) . 
which were dominated by high-mode neutrino-driven convec- 
tion, we find that the low-mode unstable models presented here 
show significantly (i.e. up to 40%) higher initial metal clump 
velocities. These high velocities, in turn, keep the timescale of 
clump propagation through the He core short - in particular 
shorter than the timescale of reverse shock formation near the 
He/H interface. Hence the fastest clumps in the low-mode mod- 
els never have to interact with this reverse shock, and are thus 
not slowed down, in marked contrast to the high-mode case 
(lKifonidisetalJ2003h. 

The occurrence of a strong Richtmyer-Meshkov instability 
at the He/H interface is of crucial importance for producing the 
strong inward mixing of hydrogen and the global anisotropy 
that is present during the late expansion of the ejecta. The early 
deposition of vorticity into the interface layer by the initially 
aspherical shock (at t ^ 100 s) is responsible for the fact that 
a late-time global anisotropy of the ejecta can develop (at t ^ 
10000 s), although the shock itself has long become spherical 
by that time. 

We find that the large-scale anisotropy imprinted on the 
ejecta and the early shock by low-mode hydrodynamic in- 
stabilities during the shock revival phase, and the ensuing 
Richtmyer-Meshkov instability at the He/H interface are all 
that is needed to reproduce the major observational features 
of SN 1987 A. In other words, mechanisms like jet formation, 
anisotropic neutrino emission, and rotation of the progenitor 
are not required to explain the observational characteristics we 
have summarized in Sect. 13 provided that due attention is paid 
to account for all hydrodynamic instabilities that occur in the 
problem and for their mutual interaction. 

A coherent picture of supernova explosions, powered by 
the neutrino-heating mechanism, thus seems to emerge, in 
which non-radial instabilities of the accretion shock like the 
advective - acous tic instability of lFoglizzo & Tagged (l2000l) and 
^Foglizzo* ('2002^ (se e in particular 'poglizz o & Galettil l2003[ 
kFoglizzo et al. 2005t lOhnishi et al. 2005), or the stalled ac- 
cretion shock instability ^(SAS I) of lBlondin etalJ (l2003h and 
iBlo ndin & Mezzacappa (2005"), may be of significant impor- 
tance for globally distorting the supernova shock during the 
shock revival phase and for accelerating th e neutron star to ve- 
locitie s in excess of several hundred km/s ('Sche ck et al.ll2004l 
l2006h . Here we have demonstrated that such an early non- 



radial, low-mode shock instability may also be of crucial im- 
portance for the observational appearance of a supernova, be- 
cause it can provide the shock anisotropy which is responsible 
for the very strong Richtmyer-Meshkov instability that we ob- 
serve at the He/H interface. Judging from the results of our 
simulations, this latter instability also deserves much more ap- 
preciation by supernova modelers than it has received so far. 

Although our simulations make use of a polar grid and the 
assumption of axisymmetry, a fact which we critically assess in 
Sect. 13.21 the development of the discussed low-mode instabil- 
ities and the corresponding explosion asymmetries is neither a 
numerical artifact, nor a consequence of nonconvergence of the 
simulations due to a lack of resolution, or even code deficien- 
cies. The variability of the early explo sion asymmetr ies found 
in the present study, and in the work of lScheck et al ( 2006), is 
the result of an extremely nonlinear, chaotic growth of the ini- 
tially imposed small random seed perturbations. No numerical 
convergence of simulations (in the mathematical sense) can be 
expected in such a situation, especially for morphological as- 
pects like the shock deformation or the matter distribution in 
the unstable layer. Such aspects diff'er even between runs of ex- 
actly the same model, if these runs are carried out on diff'erent 
computers with slightly diff'erent 64 bit round-off' behaviour. 
However, integral quantities like the explosion timescale and 
energy, or the neutron star mass, show little diff'erences between 
simulations which make use of the same boundary conditions, 
in spite of large variations of the explosion geometr y. Based on 
a huge set of models this is discussed in detail in S check et al.l 
(2006). Therefore we do not expect that the findings presented 
here will be challenged by future simulations with higher grid 
resolution or diff'erent numerical schemes. 

Furthermore, the occurrence and the final extent of the mix- 
ing due to the Richtmyer-Meshkov instability is fairly robust 
against the range of early explosion asymmetries that was ex- 
plored in the present pa per. All three low-mod e models that we 
have computed for the Wooslev et al.! ('l988) progenitor (i.e.. 
Models bl8a, bl8b, and b23a) show a comparable extent for 
the inward mixing of hydrogen and the outward mixing of met- 
als in mass (for a comparison between Models bl8b and b23a, 
see Fig.|6ll. A larger sample of simulations is, however, required 
for a more reliable assessment of this issue. 

The non-rotating models that we have presented here pre- 
dict non-rotating neutron stars (unless an off'-center kick im- 
parts also a sizeable spin to the neutron star). This is in agree- 
ment with the fact that to date no pulsar has been found in 
SN 1987 A. However, this point could be satisfied equally well 
by a rotating neutron star with a weak surface magnetic field. 
Such an object would not produce significant pulsar activity, 
and its thermal cooling signature falls below the present detec- 
tion limit. It remains to be seen in future simulations, though, 
whether rotating models can be constructed that are also able 
to reproduce the other observational items that we have sum- 
marized above. 

A pulsar would, of course, also be lacking if a black hole 
had formed in SN 1987 A. This option, however, is quite un- 
likely for an 18-20 M© progenitor and with the core struc- 
ture of current progenitor models would require an unaccept- 
ably soft nuclear equation of state, which seems to be excluded 
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by the recent detection of a neutron star with a well-measured 
mass of 2.1 ± 0.2 Mq (Nice et al. 2005). 

While we are well aware of the fact that comparably well 
resolved three-dimensional models and detailed lightcurves, 
spectra, and polarization data still remain to be calculated, we 
think that already the present results provide support to the hy- 
pothesis that SN 1987 A was an explosion which fits into the 
framework of the neutrino-driven explosion paradigm, and that 
its understanding does not require the assumption of (more) 
controversial physics. 

In the future we plan to run long-time simulations for Type 
Ib,c explosions, in which the He/H interface is absent and the 
large anisotropy of the metal and/or helium core which we 
found in our models will be exposed to the observer from the 
outset. We believe that the global asymmetries imprinted on the 
ejecta by low-mode hydrodynamic instabilities in the neutrino- 
heated postshock layer may also off'er an explanation of the 
large-scale anisotropy of the element distribution visible in the 
Cassiopeia A supernova remnant ( Hwang et aL.2004 ). Our re- 
sults in Figs.[8land|9l and the results of Paper I, suggest that in 
this case the velocities of iron-group material can be expected 
to be even higher for a given energy of the explosion. 

Appendix A: The AUSM+ flux 

We consider the construction of a numerical flux function for 
the system of Euler equations 

dU dF{U) 



dt 



dx 



= 0, 
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p, w, and p have their usual meanings. 
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is the (specific) total energy, e the internal energy, and H the 
total enthalpy. We shall only summarize the calculation of the 
first order accurate AUSM+ flux for ID Cartesian geometry 
here. The extension of the method to higher orders of accuracy, 
curvilinear grids, and more than one spatial dimension is done 
usi ng standard techniques, (see e.g. .Colella & Woodward 19841 
and lLioul 19961) and will not b e explained i n detail here. 

The AUSM+ scheme of ^Lioul (ll£96^ is a sequel to the 
Advection Upstream Splitting Me thod (AUSM) that was in- 
troduced bv lLiou & Steff'enI (1 19931) . A central idea in these sec- 
ond generation flux vector split schemes is that the advective 
and acoustic waves of the system (lA.U describe two physically 
distinct processes and should be treated as such by the numer- 
ical scheme. To achieve this, the following ansatz is made in 
AUSM+ to calculate the numerical flux, /j+1/2, at the interface 
J + 1/2 between zones j and 7 + 1 



/7+1/2 = mj+ii2aj+ii2<5>j+ii2 + 



^ 

Pj+m 
j 



(A.4) 



i.e. /j+1/2 is written as the sum of a convective and a pres- 
sure term. In Eq. (IA.4ft the quantities my+1/2, aj+1/2, Pj+1/2 and 
Oy+1/2 are interface values of the Mach number, sound speed, 
pressure and the modified state vector O = (p,pu,pH)^ , which 
contains the total enthalpy, //, instead of the total energy, E. 

The scheme proceeds by accounting for the eff'ects of right- 
going and leftgoing nonlinear waves (with eigenvalues u ± a 
of the flux Jacobian dF(U)/dU) by a proper definition of the 
quantities m^+1/2 and Pj+i/2- The eff'ects of the linearly degen- 
erate field with eigenvalue u are finally taken care of by a sim- 
ple upwind selection of ^j+i/2- The choice of aj+1/2 mainly 
aff'ects the resolution of shocks (and is discussed in more de- 
tail belo w). Th e formulae for m^+1/2, Pj+1/2, ^7+1/2, and 0^+1/2 
given guarantee up winding, an exact resolution 

of stationary shocks and of both stationary and moving con- 
tact discontinuities. Furthermore they ensure the positivity of 
the scheme, and most important for the present application, the 
avoidance of the "odd-even-decoupling" and carbuncle phe- 
nomena in multidimensional applications (see Liou.2000L for 
details). 

Liou's formula for the interface Mach number is 



mj^i/2 = M\Mj)^M-(Mj^i), 



(A.5) 



where the contributions of the rightgoing and leftgoing waves 
have been denoted with superscripts "+" and respectively. 
Here 

Mj = Uj/aj+i/2, Mj+i = Uj+i/aj+i/2 (A.6) 

are left and right Mach numbers that are computed using the 
common interface sound speed aj+1/2. The split Mach numbers 
M- are given by the functions 



j|(M±|M|), if|M|>l, 
M-(M) = ^ 2 . 

otherwise. 



[At|(M), 
where the Ma are defined as 



(A.7) 



^t|(M) = ±i(M±l)2±yS(M2-l)^ p= 1/8 (A.8) 

(and w e have corr ected a misprint in the corresponding equa- 
tion of lLioulfr996l) . Similarly, the interface pressure is defined 



as 



Pj+i/2 = P^(Mj)pj + p-(Mj+i)pj+u 
where the split pressures are given by 

^ (1 ± sign(M)) , if|M| > 1, 



P-(M) 



otherwise. 



(A.9) 



(A. 10) 



with 



1 2 
^±(M) = -{M± \f (2 + M) ± aM^lVp- - l) , a = 3/16. 

The final upwind selection of 0^+1/2 can be written con- 
cisely in terms of the quantities 

^~j+ll2 = \ (^7+1/2 ± h7+l/2|) . (A. 11) 
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With their help the numerical flux at interface 7 + 1/2 can be 
expressed as 



/7+1/2 = aj+ii2\m)_ 



7+1/2 



[p] 




(P] 






' ^ 


pu 




pu 




> + 






j 








[ j 



What remains to be specified is the interface sound speed 
aj+i/2. In principle the simple arithmetic average 

cij+i/2 = ^ {aj + aj+i) (A.12) 

may be taken a t the e xpense of slightly smearing shock waves. 
The formula of lLioul 11994) 



with 



aj+i/2 = min(ay,ay+i), 



a = a^/ max (a*, \u\). 



(A.13) 



(A. 14) 



and the critical speed of sound, a^, calculated via the isoener- 
getic condition for an ideal gas 



H = 



1 



T^2"" 



(A.15) 



r-i z 2(r-i) 

gives sharper shocks. In particular it allows for exact resolution 
of stationary shocks (the generalization of (IA.15II to more gen- 
eral equations of state, in which the ratio of specific heats, 7, 
is not a constant, is straightforward). Moving shocks are gener- 
ally captured within one to two zones. 

Numerical experiments confirm that the scheme is substan- 
tially less viscous than older flux vector splittings. In fact we 
have found the AUSM+ flux to be even less viscous than the 
Godunov flux (obtained from the exact solution of the Riemann 
problem), which makes the former a somewhat less robust 
building block than the latter if used in higher order schemes to 
compute problems with complicated sho ck interactions (as e.g. 
thewell-known blastwaves problem of IWoodward & Colellai 
Il984h . For the present application, however, in which we em- 
ployed the AUSM+ flux only near a single strong shock (where 
furthermore the flattening algorithm of PPM reduced the spa- 
tial reconstruction to first order accuracy), the method worked 
well. 
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